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In  the  main  program  listing: 

page  72     Change  card  number  155  to  read 

DO  601  N=l,  NUMNP  155 

page  73    Remove  card  preceding  card  number  180  that  reads 
189   CALL  SYMSOL  (1) 

page  74   Change  statement  number  on  card  number  203  to  read 

1020  FORMAT  (6I5,D10.3)  203 

In  SUBROUTINE  FORM: 

page  76   Insert  a  missing  card  between  card  numbers  285 
and  287   to  read 

CV  m   SPHT(MTYPE)  *  CFUNC(MTYPE,TMEAN)  286 

In  the  MAIN  Program  listing: 

page  72   Change  Card  Number  138  to  read 

TEMP=TEMPB(N)  -  T<2>  138 

Change  the  +  signs  to  -  signs  on  the 
following  Card  Numbers: 


144  ~l(Note:   +  Q   is  OUT  ) 

145  J 


146 
147 
149 
151 
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I.   INTRODUCTION 

In  this  chapter  the  reasons  and  advantages  of  the  finite  element 
analysis  method  are  given.   The  purpose  of  this  presentation  and  the  ob- 
jectives of  the  study  are  discussed. 
A.   REASONS  FOR  THE  FINITE  ELEMENT  ANALYSIS 

With  the  ever  increasing  complexity  of  problems  occuring  from  engi- 
neering situations  in  design  and  analysis,  there  is  a  rising  need  for  an 
approximate  method  of  solution  for  these  problems.   This  technique  must 
be  flexible  enough  to  encompass  a  large  variety  of  problems  and  still  give 
accurate  answers  to  the  simplest  of  situations.   Since  most  approximate 
methods  of  analysis  of  complex  engineering  problems  involve  many  tedious 
calculations,  it  would  enhance  the  usefulness  of  the  technique  if  it  were 
capable  of  being  programed  for  computer  solution.   This  would  relieve 
the  engineer  of  the  tedious  burden  of  hand  calculated  solutions  and  leave 
him  more  time  to  devote  to  the  job  of  engineering. 

The  finite  element  method  of  analysis  is  an  approximate  method  of 
analyzing  engineering  problems  which  meet  these  specifications.   It  has 
a  high  degree  of  flexibility  and  can  be  readily  applied  to  computer  sol- 
utions to  obtain  rapid  solutions  to  complex  problems. 
B0   OBJECTIVES  OF  THIS  STUDY 

The  major  objective  of  this  study  is  a  computer  program  using  the 
finite  element  method  of  analysis,  for  the  solution  of  two-dimensional 
unsteady  heat  conduction  problems.   In  fulfilling  this  goal  there  were 
several  minor  objectives  accomplished  along  the  way.   The  original  pro- 
gram was  written  for  the  IBM  7094  system.   It  was  converted  for  operation 
on  the  IBM/OS  360  computer  at  the  Naval  Postgraduate  School  Computer  Center, 
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Involved  in  this  conversion  was  the  changing  of  control  cards  and  a  few 
FORTRAN  IV  language  differences.   Also  faster  random  access  disks  were 
substituted  for  the  magnetic  tapes  originally  used  for  temporary  storage. 
The  input  and  output  formats  were  overhauled  to  make  better  use  of  space 
on  the  printed  output  page.  More  extensive  labels  were  added  to  the 
printed  output  in  order  to  ease  the  use  of  the  program  by  persons  not 
familiar  with  the  finite  element  technique. 
C.   PURPOSE  OF  THE  PRESENTATION 

The  purpose  of  this  presentation  is  to  make  available  a  computer 
program  using  the  finite  element  method  of  analysis  for  the  solution  of 
unsteady  heat  conduction  problems.   In  doing  so  an  introduction  into  the 
formulation  of  the  finite  element  method  is  given  for  those  not  familiar 
with  the  technique,  to  give  an  indication  of  the  power  of  this  type  of 
analysis.   A  description  of  the  components  of  the  program  is  included  to 
enable  the  user  to  better  understand  the  logic  of  the  computer  solution 
and  increase  the  user's  ability  to  formulate  problems  to  which  the  pro- 
gram can  be  applied.   Specific  instructions  on  formulation  of  problems 
and  loading  of  input  data  is  made  available  in  a  chapter  on  user  infor- 
mation.  Finally,  as  is  customary  in  works  of  this  type,  recommendations 
for  future  modifications  of  the  program  and  possible  applications  are 
inc  lud  ed , 
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II.   THE  FINITE  ELEMENT  FORMULATION 

In  this  chapter  the  finite  element  method  is  formulated.   The  idea 
of  a  functional  of  a  function  and  the  variation  of  a  functional  are  intro- 
duced.  These  principles  are  developed  into  an  approximate  solution  of 
the  two-dimensional  unsteady  heat  conduction  equation  using  a  matrix 
iterative  solution. 

A.   FORMULATION  OF  THE  VARIATIONAL  METHOD 

If  a  functional  of  a  temperature  field  can  be  found  such  that  its 
extremum  yields  an  equation  describing  unsteady  heat  conduction, 

and  an  equation  describing  the  heat  flux  across  the  surface  of  the  region, 

ft*  fcjT>j  =  *iy  (2-2) 

then  the  variation  of  the  functional  will  yield  a  set  of  first  order, 
ordinary  differential  equations  in  terms  of  the  nodal  temperatures.   Using 
the  method  of  tensor  notation,  a  comma  denotes  differentiation  with  res- 
pect to  the  following  subscript  and  repeated  subscripts  imply  summation. 

Let  TT  be  a  functional  of  the  temperature  field  T(x,y,z,t)  and  the 
first  time  derivative  of  the  temperature  field  T(x,y,z,t),  be  defined  by 

TTCT,  t)  -taffc  kq  ^fcT  t-  iT)dV  +  J"  hfljas  (2.3) 

The  variation  of   fl (T>T)  with  respect  to  T  (with  T  held  constant) 
is  given  by 

iTT=iTr(,T+6/l,T)  (2.4) 
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where  £.   is  a  small  parameter  and   A  is  one  of  a  family  of  functions. 
A  is  zero  on  the  boundary  of  the  region  and  arbitrary  elsewhere.   An 
extremum  of  the  functional  "TT(T,T)  implies  that   «TT(T,T)  equals  zero,  i.e., 


-  o 


e  =o 


First 


TI  0>  tX,t)  =    (i(r+frA),.  k  tj  (T+  fe  AU 
+  ec(J<-frX)T-t(T^eX))ay-UiV(T+6^dS 


Using  the  identity 


if  the  conductivity  tensor  is  symmetric  i.e.  kij  «=  kji.   Then 

3n(r+*X,t)  _  K|r+*A),tllyX|    "  k-t](T4fcX),..X 

The  volume  integral  H  Nj  RtjA/.»OY  can  be  transformed  int 


(2-5) 


(2-6) 


(2-7) 


o  a  surface 


integral. 


(2-8) 
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When  (2-7)  is  evaluated  at  £   «=  0 


or 


The  extremum  of  Tl       f  ^TT-  0  )reduces  to 


(2-9) 


kcjT?jC  -  eci  -%  (2-10) 

in  the  region  V.   Also  on  the  surface  S 

^kyXj  *  n-^  (2-ii) 

These  two  equations  (2-10,  2-11)  are  the  unsteady  heat  conduction  equation 
and  the  boundary  flux  equation.   This  verifies  that  a  function  T  which 
yields  an  extremum  of  the  functional  defined  by  (2-3)  satisfies  both  the 
unsteady  heat  conduction  equation  in  the  region  V  and  the  boundary  flux 
equation  on  the  surface  S    Qjj  . 

B.   FORMULATION  OF  THE  FINITE  ELEMENT  PROBLEM 

If  the  choice  of  T  and  T  is  restricted  so  that  they  are  represented 
by  certain  constants  which  are  obtained  in  their  formulation,  the  func- 
tional IT  becomes  a  real-valued  function.   In  this  case  \\    shall  be 
^QlFjTI'i/   where  {T}   an°l  iTj  are  vectors  of  nodal  point  values 
of  (t|  and  yTi£%  The  extremum  of  TT^VTv  ,  Vl"}  )  must  now  be  found  which 


Numbers  in  brackets  refer  to  similarly  numbered  items  in  Bibli- 
ography. 
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implies 

■frdTUtfO      =     O  (2-12) 

From  this  point  on  the  finite  element  formulation  will  be  considered 
in  two  dimensions  only.   The  region  V  will  be  divided  into  a  number  plane 
triangular  elements.   The  functions  T  and  T  can  be  uniquely  described 
within  each  of  these  elements  by  the  values  of  the  function  at  each  of 
the  nodal  points  of  the  element  T^,  Tj  ,  Tm,  and  linear  function  of  the 
coprdinates  of  the  element. 

Let  the  temperature  field  in  an  element  be  given  by 

TUv,-t-)  =  <95cx,v)>|A]^}     ?=sj,m  (2.13) 

and  the  time  rate  of  temperature  change  given  by 

TU.v,t)=<tf(xv)>[A](t?U)]  (2.14) 

\X>y    is  a  vector  which  specifies  the  spacial  approximation  of  T  and  T. 
The  matrix  [Aj  is  a  constant  and  is  defined  by  the  above  relationships. 
Its  exact  value  and  derivation  will  be  discussed  in  detail  in  section  II, 
E. 

The  temperature  gradient  T,*  ( ft  n  x,y)  can  be  expressed  in  terms  of 
the  nodal  temperatures  by 

In  the  following  development  the  subscript  f    ,  denoting  nodal  point 
values,  will  be  dropped  and  will  be  assumed. 
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Writing  the  functional  TT  in  terms  of  nodal  point  values, 


T,,vT.  \  P.    X  JT 


-{T}  [A]<^>i)dv    -]iTjW4ftWitIJS  (2-16) 

Taking  the  variation  of  TT  with  respect  to  \t}   and  setting  it  equal  to 


zero 


*ttUtu+})=[«](ti  -grh&ift-fts-o    (2.17) 


where 


T..,.l 


[*]«      WOT  W6*3WJv  (2.18) 


Jy 

[tl  =  J  iWVd 


(2-19) 


V  (2-20) 

fori  - 1  wWwm<» 

C.   BOUNDARY  CONDITIONS 

There  are  four  types  of  boundary  conditions  to  be  considered.   (1) 

specified  temperature,  (2)  specified  flux,  (3)  convection  boundary  layer, 

and  (4)  adiabatic.   For  a  specified  temperature  boundary  condition  T^  is 

constant  over  boundary  segment  S-^.   For  a  specified  flux  q  t=  q  over 

boundary  segment  S2.   For  a  convection  boundary  layer  q  a  h(Tf  -  T0)  over 

boundary  segment  Sg.   Finally  for  an  adiabatic  boundary  condition  q  ■  0 

over  boundary  segment  S4. 
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where 


The  boundary  integral  (2-21)  becomes 

{  V}    '-      [AJ<0>\}%]  dg  (2-23) 


[M]  -U  [a]V  <#>[*]  <te  (2.2« 


% 


5h*}  -"   kX   [AlWdS  (2-25) 

The  integral  on  segment  S  is  zero  because  the  variation  of  the  functional 
was  specified  as  zero  on  the  boundary.   The  integral  on  segment  S^  is 
zero  since  there  is  no  heat  flow  across  the  boundary. 

To  obtain  the  extremum  of  the  functional  TT  the  nodal  point  tempera- 
tures  yl]  which  satisfy  the  following  matrix  equation  must  be  found   |~6j . 

(W  -U)\t]  +  fc](t]  -  {%*}  -$>  [ir]  --  o  (2.26) 

D.      THE  MATRIX  EQUATION 

Let  the  time  variable  be  ti  such  that  t^  =  i  A"t     i  =  0,1,2,... 
For  this  development   i  will  refer  to  time  increments. 

(2-27) 


M  -  [K]-tH] 


Equation  (2-17)  can  be  written  as 


(2-28) 
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Assuming  that  \  T  r  is  constant  over  the  interval  t^^  t<j£.tj  ,  ,,  i.e., 

mi-  cm,41  -iykv^ 

1 T j  ,       can  be  obtained  by  a  Taylor's  expansion  around   t  ■  t£. 

•     iTlfc,    «  iTh  *  ^U  *    f  tffy  (2-30) 

m . ,  MTl;  ► *tif i :  *  fttlii-iili)         (2-3D 
ItW  =  {r]t  ^(irL,  -  It}. )  (2-32) 

Using  equations  (2-29)  and  (2-32)   ^T].^,  and  {Tl<        can  be  determined 
in  terms  of  (T^,   and    (l],  .   Solving  (2-32)  for  {t}.    gives 


ftL,"  Balp3u-mi)-it}t 


(2-33) 


Substituting  this  into  (2-28)  gives 


K  -f&Ea  M*.  ■  ffta  +  E3(AIT}C  t(T]L)     ^3*) 


Also  from  equation   (2-28) 
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Substituting  (2-34)  into  (2-35)  gives 


([k]4J3)H^H%K  4fk  +tfl 


(2-36) 


or 


([k]  4i|ty^) = mirh+mA^ 


(2-37) 


In  a  multi-element  body  the  element  matrix  equations  must  be  assembled 

into  a  single  matrix  equation.   The  assembly  and  solution  of  the  equation 

is  a  complex  task  but  it  can  be  greatly  simplified  by  the  use  of  a  com- 
puter [e]. 

E.   FORMULATION  OF  THE  ELEMENT  MATRIX 

If  the  region  to  be  analyzed  is  divided  into  plane  triangular  elements 
then  a  typical  element  would  appear  like  the  following  figure. 


Figure  2-1  Triangular  Element 
The  temperature  in  an  element  can  be  described  in  terms  of  a  linear  com- 
bination of  the  coordinates  of  the  elements  nodal  points  and  the  values 
of  the  temperature  at  each  of  the  nodal  points. 

Assuming  that  the  temperature  distribution  is  of  the  form 


"nx,Y,o  =  o(,a)  -v^cox  *a<3t-oy 


(2-38) 
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The  temperature  at  each  of  the  nodal  points  can  be  written  as 


Tj  "  °*>      *         ^j^2-   f  Yj0^ 


(2-39) 


Tm   -    «,   +     Xm*t    V   ym<^ 


Solving   for    the  coefficients    °^ -i ,   °Co»    an^        °Vo   using   Cramer's   rule  the 
following   results   are  obtained. 


Lglt 


Note: 


I  X;  HI! 


=    £A 


is  the  area  of  the  triangular  element. 


(2-40) 


<*.  = 


Tt 

*i 

yj 

1* 

*j 

i!i 

"B. 

x* 

1m 

«t  = 


; 

~:   * 

/ 

1^ 

1 

T     V 

2.A 


IA 


"fe- 


ZA 
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Recalling  that  the  temperature  distribution  was  previously  described  by 
the  expression 


I 


=  <^>[Al[T] 


(2-42) 


)t\      is  the  vector  of  the  nodal  temperatures  and  JAj  is  the  spacial 
dependence  matrix  which  is  a  constant.   Let 


<0>   =    <|      X     v> 


(2-43) 


hi- 


then 


[A]- 


<n* 


°0'       °<3^    ^33 


=  <0>[A]{T} 
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so 


r 


*c 


y 


r-6 


Therefore,  comparing  (2-43a)  with  (2-41) 


l/fc. 


°*»3 


2J 


x  . 

J 

Xfrvy 

V 

X  W\ 

7, 

1: 

X' 
J 

I         7; 


^V^-U7j 


^**  /c  Al  I  w 


*>i-*i)t 


h  ~ 


V 


m 


°Vzz      — 


7* 


'iv»  ~  '<• 


°*Z3     - 


°<3< 


$ 


Xj 


** 


Y;  -  // 


/C  yyx      -     *j' 
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li-  X 


wv 


'S3 


A: 


Xj  -1c 


The  spacial  dependence  matrix  is  then  defined  by 


[A]     = 


ZA 


y,  -v      y .  v,       v._v  / 


A  ^A    *j        Xv 


Xj-v: 


(2-44) 


DO. 


F.   THE  METHOD  OF  LUMPING  ELEMENT  PROPERTIES 

When  working  with  a  computer  it  is  usually  more  convenient  to  work 
with  a  quadrilateral  element.   The  quadrilateral  is  divided  into  four 
triangles  with  the  center  nodal  point  common  to  each  of  the  triangles, 
(see  figure  2-2). 


->  X 


Figure  2-2   Quadrilateral  Element 
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The  matrix  equation  (2-28)  can  be  assembled  for  this  element  in  the 
following  form 


Kll  K12  K13  K14  K15 


K21  K22  K23  K24  K25 

K31  K32  K33  K34  K35 

K41  K42  K43  K44  K45 
K51  K52  K53  K54  K55 
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Partitioning  this  equation  results  in 


^5j>  K55 


H  + 


(C5j>  C55 


H 


(2-45) 


i,j  B  1,2,3,4 

This  equation  can  be  multiplied  out  to  yield  two  equations 

h^m  +  [K«]  Ts  +  [cij]l*A  +  M  fs  ■  £M  (2-46) 

<K5j/lT4     +      K55T5       +         <C5j>l?i]       +       C55?5       ■       f5  »"*» 

The  specific   heat  matrix     ]  C  1       is  approximated  by    lumping   the  heat   ca- 
pacities  at   the   four   external  nodes.      Thus     |Cj    becomes   a  diagonal  matrix 
(4x4)    and   C55   =  0.      Therefore    (2-46)    can  be  written  as 

[Kij]   lTi]        +         fKi5]        T5  +         [C±I]{H        ^  {fi}  <2-48> 


and    (2-47)    as 


<K5j)    {Ti]        +         K55T5 


(2-49) 


Solving    (2-49)    for  T,.    and   substituting  into    (2-48)    results   in 
LKii]    {  Ti]     +    (Ki5]  K55  (  f5   "  <K5j>[Ti"}  )    +     [Cij]  {  fi]     =  (f i]        ^2-5°) 
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I 


or 


(  &$    "   Kl  K55  (\)  >   {Ti}     +     CCij](*i]    "  (f i]   "   {Ki5}  K55f5      <2"51) 

Equation    (2-51)    is   analogous   to   equation    (2-28)    except   that    JKJ     and      fc"l 
are  now   (4x4)   matrices   and  Sf\     is   a    (4x4)   vector      [6j   . 
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III.   PROGRAM  STRUCTURE 

This  chapter  describes  the  general  composition  of  the  program  and 
gives  a  description  of  the  function  of  each  section  of  the  program. 

A.  PROGRAM  IDENTIFICATION 

The  program  was  written  by  Robert  E.  Nickell  while  in  the  Department  of 
Civil  Engineering,  University  of  California  at  Berkeley.   It  was  formerly 
identified  as  W2498,W035--2--57343,  Nickell.   It  was  written  in  July  1968 
and  was  revised  by  Allan  B.  Chaloupka  in  May  1969.   Formulation  was 
based  upon  a  variational  principle  derived  by  M.  E.  Gurtin  of  Brown 
University  [Jt J . 

The  program  has  been  renamed  DFETHC-Finite  Element  Transient  Heat 
Conduction. 

B.  PURPOSE 

The  purpose  of  the  program  is  to  produce  high  speed  solutions  to 
transient  heat  conduction  problems.   These  problems  may  involve  plane  or 
axisymmetric  geometry,  non-homogeneous  anisotropic  material  configura- 
tions, temperature  dependent  material  properties,  and  time  dependent 
boundary  conditions. 

C.  PROGRAMING  INFORMATION 

The  program  was  originally  written  in  FORTRAN  IV  for  the  IBM  7094 
computer.   It  was  converted  for  use  on  the  IBM  OS/360  Model  67  computer. 
The  program  was  also  rewritten  in  double  precision  (REAL  *8) . 

D.  PROGRAM  CAPACITY 

The  size  problem  that  can  be  analyzed  by  this  program  is  subject  to 
the  following  limitations" 

Maximum  number  of  nodal  points  500 
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Maximum  number  of  quadrilateral  elements 490 

Maximum  number  of  materials 25 

Maximum  difference  of  nodal  point  numbers 

for  the  same  e  lement  „  „ , 31 

These  are  not  stringent  limitations  and  they  may  be  adjusted  some- 
what by  the  user  to  fit  individual  problems.   The  limits  are  changed  by 
changing  the  appropriate  DIMENSION  and  COMMON  statements  and  one  card  in 
the  segment  that  calculates  the  half -band  width.   The  only  real  limita- 
tion placed  upon  the  size  of  the  problem  that  can  be  analyzed  is  the 
amount  of  core  storage  available  in  the  computer  used. 
E.   GENERAL  STRUCTURE 

The  program  consists  of  a  main  program  and  several  subprograms  used 
for  repetitive  operations.   The  main  program  can  be  divided  into  three 
parts.   The  first  part  is  concerned  with  reading  into  the  program  and 
printing  out  the  data  describing  each  problem.   It  is  portrayed  by  the 
flow  diagram  in  figure  3-1. 

The  second  portion  of  the  main  program  is  concerned  with  the  routing 
of  operations  according  to  certain  logical  variables.   These  variables 
describe  four  possible  combinations.   They  are  constant  material  proper- 
ties and  constant  boundary  conditions,  temperature  dependent  material 
properties  and  constant  boundary  conditions,  constant  material  proper- 
ties and  time  dependent  boundary  conditions,  and  temperature  dependent 
material  properties  and  time  dependent  boundary  conditions.   Each  of  these 
possibilities  is  portrayed  by  a  flow  diagram  in  figures  3-2  through  3-5. 

The  third  portion  of  the  main  program  is  concerned  with  calculating 
the  effective  forcing  vector,  solving  the  matrix  equation  and  printing 
the  temperatures  for  each  time  increment.   The  operation  is  then  directed 
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Control  Parameters 


Read  and  Print 
Material  Properties 


Read  and  Print 
Time  Sequencing  Information 


Read  or  Generate  and  Print 
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Figure  3-1  Main  Program  Flow  Diagram,  Part  I 
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Figure  3-2  Main  Program  Flow  Diagram,  Part  II 
Constant  Properties  and  Boundary  Conditions 
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Figure  3-3  Main  Program  Flow  Diagram,  Part  II 

Temperature  Dependent  Properties  and 

Constant  Boundary  Conditions 
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Figure  3^4  Main  Program  Flow  Diagram,  Part  II 
Constant  Properties  and  Time  Dependent  Boundary  Conditions 
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Figure  3-5  Main  Program  Flow  Diagram,  Part  II 

Temperature  Dependent  Properties  and 

Time  Dependent  Boundary  Conditions 
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back  to  the  beginning  of  the  second  section  of  the  main  program  to  start 
calculations  for  the  next  time  increment.   The  third  portion  of  the  main 
program  is  portrayed  by  a  flow  diagram  in  figure  3-6. 
F.   MAIN  PROGRAM 

The  main  program  is  composed  of  a  number  of  individual  steps  which 
comprise  the  finite  element  solution  technique,  and  a  system  for  the  in- 
put and  output  of  problem  information.   Not  included  in  the  main  program 
are  the  matrix  assembly  processes,  the  matrix  equation  solving  processes, 
and  the  boundary  condition  input  routines.   These  are  included  in  sub- 
programs FORM,SYMSOL,TFUNC,CFUNC,FCBC,FTBC,  and  FFBC. 
1.   Control  Information 

The  first  information  to  be  read  into  the  program  is  the  control 
information.   This  information  sets  the  basic  variables  for  each  of  the 
steps  of  the  program  and  also  routes  the  problem  to  various  sections  of 
the  program.   This  information  includes  the  number  of  nodal  points,  the 
number  of  elements,  the  number  of  convection  boundary  conditions,  the  num- 
ber of  materials,  the  number  of  time  sequences,  whether  the  problem  is  in 
plane  or  axisymmetric  geometry,  the  output  print  interval,  the  system  of 
units,  whether  a  final  spacial  temperature  distribution  is  desired,  and 
the  reference  temperature  of  the  body. 

Also  read  in  with  the  control  information  are  two  titles.   The 
first  is  a  72  space  title  that  can  be  used  by  the  user  to  label  each 
problem.   The  second  is  a  48  space  title  that  is  used  to  input  the  labels 
used  with  the  system  of  units  chosen. 

The  appropriate  titles  and  control  information  are  printed  on  the 
first  page  of  the  printed  output,  for  each  problem. 
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Figure  3-6  Main  Program  Flow  Diagram,  Part  III 
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2.  Material  Properties 

The  thermal  conductivity,  density,  specific  heat,  and  heat  gener- 
ation properties  for  each  material  are  read  into  the  program.   The  materi- 
al properties  are  then  printed  in  order  on  the  printed  output  page  and 
labeled  with  the  appropriate  units.   Since  the  program  is  for  anisotropic 
materials,  kxx,kVy,  and  kxy  are  read  in  order  to  form  the  thermal  con- 
ductivity tensor.   The  specific  heat  and  density  are  read  in  as  a  product 
in  order  to  reduce  the  number  of  variables  in  the  program. 

3.  Time  Sequencing  Information 

In  this  section  the  time  sequencing  information  is  read  in  for 
each  sequencing  scheme.   This  information  consists  of  two  logical  varia- 
bles and  two  numerical  variables.   The  logical  variables  indicate  whether 
there  are  temperature  dependent  properties  and  if  there  are  time  dependent 
boundary  conditions.   The  two  remaining  variables  describe  the  number 
of  time  steps  and  the  size  of  each  time  increment.   The  time  sequencing 
information  is  printed  on  the  printed  output  page  with  appropriate  labels. 

4.  Nodal  Point  Information 

Ordinarily  this  segment  of  the  program  will  read  the  x,y  coordinates 
and  initial  temperature  of  each  nodal  point  and  print  the  information 
with  the  appropriate  dimensions.   Because  of  the  possibility  of  a  large 
number  of  nodal  points  (500)  and  the  equally  large  possibility  of  human 
error  in  punching  data  cards  there  is  a  segment  of  the  program  that  will 
generate  nodal  point  information.   This  segment  has  a  few  limitations, 
which  will  be  discussed  in  chapter  IV,  but  it  is  usually  able  to  generate 
almost  all  of  the  interior  nodal  points.   As  an  example  in  the  configura- 
tion shown  in  figure  3-7  the  nodal  points  indicated  by  squares  were  read 
into  the  program  and  the  remaining  nodal  points  were  generated  by  the 
program. 
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Figure  3-7  Nodal  Point  and  Element  Generation 
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5.  Element  Information 

The  purpose  of  this  segment  is  about  the  same  as  the  previous 
one.   The  number  of  the  element,  the  number  of  each  of  its  nodes,  the 
type  of  material,  and  the  amount  of  heat  generation  within  the  element 
are  read  into  the  program  and  printed  out  with  appropriate  labels.   Again 
this  segment  has  the  capability  to  facilitate  the  loading  of  data  by 
generating  element  information.   This  segment  also  has  some  limitations 
which  will  be  discussed  in  chapter  IV,  but  it  is  usually  able  to  generate 
most  of  the  elements  in  a  given  configuration.   As  an  example,  in  figure 
3-7  all  the  elements  were  generated  except  those  that  are  starred. 

Also  included  in  this  segment  is  a  routine  that  calculates  the 
half -band  width  of  the  problem.   This  is  an  important  control  parameter 
that  is  used  later  in  the  program  to  assemble  the  matrix  equation. 

6.  Time  Steps 

The  time  step  portion  of  the  program  includes  two  large  loops. 
The  first  loop  executes  the  program  through  each  time  sequence  of  the 
problem.   The  second  loop  executes  the  program  through  the  required  num- 
ber of  time  increments  for  each  time  sequence.   Inside  these  two  loops 
the  two  logical  time  sequencing  variables  and  the  loop  counters  control 
the  routing  of  the  problem. 

The  two  logical  variables  control  the  choice  of  variable  material 
properties  and  boundary  conditions.   If  the  material  properties  and  bound- 
ary conditions  are  constant  then  the  coefficient  matrix  and  forcing  vector 
are  calculated  only  once.   The  matrix  equation  is  then  solved  for  suc- 
cessive time  increments  using  only  the  original  matrices.   If  variable 
material  properties  are  involved  the  effective  coefficient  matrix  must 
be  reapplied  since  the  coefficient  matrix  is  reformed  from  zero  values. 
If  time  dependent  boundary  conditions  are  involved  then  the  contribution 


of  the  boundary  conditions  to  the  coefficient  matrix  must  be  reapplied. 
The  coefficient  matrix  from  the  original  material  properties  is  kept  in 
temporary  storage.   When  both  time  dependent  boundary  conditions  and 
temperature  dependent  properties  are  used  then  the  coefficient  matrix 
and  forcing  vectors  must  be  reformed  for  each  time  increment  based  upon 
the  immediate  time  and  temperature  values  and  initial  properties.   The 
data  for  temperature  dependent  properties  is  calculated  from  the  initial 
conditions  and  the  present  nodal  temperatures.   The  data  for  time  depend- 
ent boundary  conditions  is  either  read  into  the  program  with  each  time 
increment  or  generated  within  the  program. 

Different  time  sequences  can  be  run  for  the  same  problem.   This 
enables  the  user  to  speed  or  slow  the  time  change  during  a  particular 
segment  of  the  problem. 

7.   Convection  Boundary  Conditions 

A  convection  boundary  condition  occurs  when  the  problem  includes 
a  convection  boundary  layer  as  part  of  its  boundary  conditions.   One  con- 
vection boundary  condition  consists  of  two  nodal  point  numbers,  the  con- 
vective  heat  transfer  coefficient  and  the  ambient  fluid  temperature.   The 
two  nodal  points  define  the  convection  boundary  for  each  condition. 

Time  varying  convection  boundary  conditions  can  be  obtained  by 
varying  the  ambient  fluid  temperature.   The  new  fluid  temperature  may  be 
read  with  each  new  time  step  or  it  may  be  generated  by  using  the  function 
subprogram  FCBC.   FCBC  is  used  in  the  same  manner  as  the  function  sub- 
programs for  temperature  dependent  material  properties.   The  argument  of 
FCBC  is  TIME  and  the  value  of  FCBC  is  currently  set  at  one.   If  tempera- 
ture dependent  material  properties  are  used  then  the  convection  boundary 
condition  section  stores  and  reapplies  the  convection  boundary  condition 
contributions  to  the  coefficient  matrix  and  the  forcing  vector. 
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8.  Temperature  and  Flux  Boundary  Conditions 

A  temperature  or  flux  boundary  condition  occurs  when  the  problem 
specifies  a  temperature  or  flux  at  a  nodal  point.   One  temperature  or  flux 
boundary  conditions  consists  of  the  nodal  point  number  and  the  value  of 
the  temperature  or  flux  at  that  nodal  point. 

Time  varying  temperature  and  flux  boundary  conditions  can  be  ob- 
tained by  varying  the  temperature  or  flux  at  each  nodal  point.   The  new 
temperature  or  flux  may  be  read  with  each  new  time  step  or  it  may  be 
generated  by  using  the  function  subprograms  FTBC  or  FFBC  for  time  depend- 
ent temperature  or  flux  respectively.   FTBC  and  FFBC  are  used  in  the 
same  manner  as  the  function  subprograms  for  temperature  dependent  materi- 
al properties.   The  argument  of  each  is  TIME.   Their  current  values  are 
both  one.   If  temperature  dependent  properties  are  used  the  temperature 
and  flux  boundary  condition  section  stores  and  reapplies  the  boundary 
conditions  to  the  coefficient  matrix  and  the  forcing  vector. 

9.  Temperature  Output 

This  segment  of  the  program  is  concerned  with  the  print  out  of 
the  temperature  solutions  for  each  nodal  point  at  each  time  step.   The 
value  of  the  time  is  printed  and  the  nodal  point  numbers  and  correspond- 
ing temperatures  for  that  time  are  printed  in  rows  below.   If  it  is  de- 
sired, the  temperature  can  be  printed  with  spacial  data  at  the  end  of 
each  time  sequence. 

G.   SUBROUTINE  FORM 

The  subroutine  FORM  is  used  to  calculate  the  effective  element  ther- 
mal conductivity,  heat  capacity,  and  heat  generation.   The  subroutine 
performs  these  operations  for  each  of  the  elements  of  the  configuration 
and,  if  temperature  dependent  properties  are  used,  for  each  time  increment. 
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First  the  mean  x,y  coordinates  and  the  mean  temperature  of  the  element 
are  calculated  by  averaging  the  values  of  each  of  the  nodes.   Using  the 
mean  x,y  coordinates  the  quadrilateral  element  is  divided  into  four  tri- 
angles.  The  conductivity  tensor,  heat  capacity,  and  heat  generation 
vectors  are  formed  for  each  of  the  four  triangular  sections.   As  these 
matrices  are  formed  they  are  assembled  into  the  larger  (5)  nodal  point 
equation  for  the  quadrilateral  element.   The  heat  capacity  and  heat  gen- 
eration vectors  are  assembled  so  that  the  fifth  term  of  each  vector  is 
dropped  when  forming  the  lumped  properties  matrix  equation.   The  center 
nodal  point  is  eliminated  from  the  thermal  conductivity  tensor  (5x5), 
according  to  equation  (2-52),  to  give  the  thermal  conductivity  tensor 
(4x4)  for  the  lumped  properties  matrix  equation.   The  thermal  conductivity 
tensor  and  heat  capacity  and  heat  generation  vectors  for  each  element 
are  returned  to  the  main  program  where  they  are  loaded  into  the  coeffi- 
cient matrix  and  vectors  used  in  the  solution  of  the  total  matrix  equa- 
tion for  each  time  step.   For  problems  where  the  properties  of  the 
materials  are  not  temperature  dependent  the  thermal  conductivity  tensor 
and  the  heat  capacity  and  heat  generation  vectors  are  calculated  only  at 
the  beginning  of  the  problem.   If  temperature  dependent  material  proper- 
ties are  included  these  matrices  are  calculated  at  the  beginning  of  each 
time  step. 

H.   SUBROUTINE  SYMSOL 

The  subroutine  SYMSOL  in  conjunction  with  the  loading  methods  of  the 
main  program  is  used  to  solve  the  finite  element  matrix  equation.   It 
consists  of  two  parts.   The  first  part  reduces  the  coefficient  matrix 
to  its  final  form.   The  second  part  forms  the  effective  forcing  vector 
and  solves  for  the  temperatures  at  each  of  the  nodal  points. 
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In  an  ordinary  formulation  the  final  matrix  equation  for  the  finite 
element  method  would  look  like  figure  3-8.   The  matrix  of  coefficients 
\R\    is  a  banded  symmetric  matrix.   In  the  computer,  in  order  to  conserve 
core  storage,  the  matrix  of  coefficients  is  stored  in  columns.   The  first 
column  is  the  diagonal  of  the  matrix  of  coefficients.   (See  figure  3-9). 
All  the  operations  performed  on  the  matrix  of  coefficients  are  done  on 
the  upper  half  of  the  symmetric  banded  matrix  in  its  stored  form. 

After  the  coefficient  matrix  is  formed  it  is  put  into  upper  triangular 
form.   This  yields  a  matrix  equation  of  the  form  shown  below. 


r  ^        r  ^ 


H 


_L  *w  J  \- 


r 


Figure  3-10  Matrix  Equation  in  Upper  Triangular  For 
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This  is  accomplished  by  the  first  segment  of  SYMS0L„   The  second  portion 
of  SYMSOL,  using  the  effective  forcing  vector,  starts  with  the  temperature 
at  the  bottom  of  the  nodal  temperature  vector  and  back  substitutes  to 
solve  for  the  nodal  temperatures.   These  operations  are  performed  for  each 
time  increment. 
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I.   TEMPERATURE  DEPENDENT  PROPERTIES 

There  are  two  function  subprograms  that  are  used  in  conjunction  with 
temperature  dependent  properties.   They  are  TFUNC  and  CFUNC.   TFUNC  is 
used  for  the  temperature  dependent  thermal  conductivity  and  CFUNC  is  used 
for  the  temperature  dependent  product  of  the  density  and  specific  heat. 
In  the  programs  present  configuration  each  of  these  function  subprograms 
is  prepared  for  temperature  independent  properties  and  their  values  are 
therefore  equal  to  one.   If  the  thermal  conductivity  were  temperature 
dependent  in  the  form  of:  k  =  kQ(l  +  bT)  ,  where  b  is  a  constant  coeffi- 
cient and  k0  is  the  value  of  the  thermal  conductivity  at  T  ■  0,  then 
TFUNC  would  have  to  be  modified  by  the  user  to  include  the  factor  1  +  bT. 
The  material  properties  read  into  the  program  would  then  be  zero  tempera- 
ture properties. 
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IV.   USER  INFORMATION 

This  chapter  describes  the  manner  in  which  problems  should  be  formu- 
lated.  Also  the  procedures  used  in  inputting  data  into  the  program  are 
included. 

A.  INPUT  DATA 

Data  for  each  problem  processed  by  this  program  is  assembled  in  groups. 
The  first  group  includes  the  titles  and  all  the  control  parameters*   The 
second  group  is  the  material  properties.   The  third  group  is  the  time  se- 
quencing information.   The  fourth  group  is  nodal  point  information.   The 
fifth  group  is  element  information.   The  sixth  group  is  boundary  condi- 
tions.  After  each  data  deck  there  must  be  three  blank  cards. 

B.  CONTROL  PARAMETERS 

The  following  control  parameters  are  read  into  the  program. 

NUMNP  number  of  nodal  point 

NUMEL  number  of  elements 

NUMCBC  number  of  convection  boundary  conditions 

NUMMAT  number  of  materials 

NSEQ  number  of  time  sequences 

KAT  type  of  geometry,  0  implies  plane  geometry 

INTER  output  print  interval 

UNIT  units  system 

KPUNCH  spacial  temperature  output 

JUMP  generated  temperature  of  flux  boundary  conditions 

KEY  generated  convection  boundary  conditions 

TO  initial  temperature 
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Also  included  with  the  control  data  are  two  labels.   TITLE  (18)  is  an 
arbitrary  label  that  may  be  used  to  identify  individual  problems.   It  is 
72  spaces  long.   The  second  label  is  UN(12).   UN(12)  is  an  array  of  four 
unit  labels.   Each  has  12  spaces.   The  units  must  be  printed  in  the  fol- 
lowing order:   length,  mass,  time,  and  temperature.   These  labels  must  be 
left  justified. 

C.  MATERIAL  PROPERTIES 

The  following  material  properties  are  read  into  the  program. 

XCOND       k 

xx 

YCOND       ky 

XYCOND      k 

xy 

SPHT        /?c 

HX  heat  generation 

Material  properties  are  read  in  with  one  material  to  a  card.   If  the 
material  is  isotropic  kxx  =  kyy  =0.   If  the  material  is  anisotropic  the 
three  values  will  not  be  the  same.   Each  material  is  given  a  number  ac- 
cording to  the  order  it  is  read  into  the  program. 

If  temperature  dependent  material  properties  are  used  both  logical 
variables  must  be  true.   Also  both  JUMP  and  KEY  must  be  some  number  other 
than  zero  and  FCBC,FTBC,  and  FFBC  must  be  one. 

D.  TIME  SEQUENCING  INFORMATION 

The  following  time  sequencing  information  is  read  into  the  program. 
LTAG        temperature  dependent  properties  (true  or  false) 
MTAG        time  dependent  boundary  conditions  (true  or  false) 
ITAG        number  of  time  steps 
TAG         time  increment 

The  time  sequencing  information  controls  the  operation  of  the  program 
through  the  use  of  the  logical  variables  LTAG  and  MTAG  but  it  can  also 
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be  used  to  slow  or  speed  time  during  the  problem.   An  example  would  be 
two  sequences  with  the  first  solving  for  the  temperature  60  times  during 
the  first  hour  of  the  problem  (i.e.,  every  minute).   The  second  time  se- 
quence would  solve  for  the  temperature  6  times  in  the  second  hour  of  prob- 
lem time  (i.e.,  every  10  minutes), 
E.   NODAL  POINT  CONSTRUCTION 

The  following  nodal  point  information  is  read  into  the  program  for 
nodal  point  data. 

N         number  of  the  nodal  point 

KODE       1  indicates  a  temperature  boundary  condition 

X         x  coordinate  of  the  nodal  point 

Y         y  coordinate  of  the  nodal  point 
/  T         initial  temperature  of  the  nodal  point 

The  data  for  each  nodal  point  is  put  on  one  card. 

Earlier  it  was  mentioned  that  the  program  could  generate  nodal  points. 
This  feature  has  a  few  limitations  that  must  be  observed  or  incorrect 
information  will  result.   A  nodal  point  whose  KODE  is  1  cannot  be  gener- 
ated.  Nodal  points  must  be  generated  along  straight  lines.   The  initial 
temperature  of  the  generated  nodal  points  must  be  the  same  and  equal  to  TO, 
the  initial  temperature.   To  generate  nodal  points  simply  leave  out  the 
data  cards  between  two  nodal  points  on  a  straight  line  and  the  program 
will  generate  the  information  for  the  nodal  points  in  between. 

In  constructing  the  mesh  of  nodal  points  for  a  problem  it  is  sometimes 
impossible  to  divide  the  area  to  be  analyzed  into  uniform  rectangular 
quadrilaterals.   At  times  it  will  not  be  desirable  to  do  so.   Consider 
figure  4-1.   In  preparing  a  nodal  mesh  for  this  configuration  the  first 
step  would  be  to  assume  that  by  symmetry  there  is  no  heat  flux  across  the 
center  lines.   This  reduces  the  figure  to  one  of  its  quadrants  as  in 
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figure  4-2.   The  size  of  the  quadrilateral  elements  would  vary  depending 
upon  how  interested  the  user  was  in  a  certain  area  of  the  figure.   In 
figure  4-2  the  area  around  the  curve  is  being  more  closely  examined  but 
the  mesh  size  must  also  be  small  to  enable  the  user  to  approximate  any 
curved  surface  using  a  series  of  straight  line  segments.   It  is  almost 
impossible  to  divide  an  area  with  a  curved  boundary  without  encountering 
a  right  triangular  element.   To  convert  this  to  a  quadrilateral  element 
place  a  nodal  point  in  the  middle  of  the  hypotenuse. 


Figure  4-3  Conversion  of  a  Right -Triangle 
to  a  Quadrilateral 


In  numbering  the  nodal  points  it  must  be  kept  in  mind  that  the  maxi- 
mum difference  in  nodal  point  numbers  in  one  element  is  31.   This  is  a 
limitation  set  by  the  half -band  width  of  the  computer  program.   It  is 
convenient  to  increase  the  numbering  of  the  nodal  points  and  the  elements 
in  the  same  pattern,  although  it  is  not  a  requirement. 
F.   ELEMENT  CONSTRUCTION 

The  following  element  information  is  read  into  the  program  for  ele- 
ment data. 

N  number  of  the  element 

IX(N,K)  K  t=  1  -  4    numbers  of  the  element's  nodal  points 

IX(N,5)  material  type 

HTGEN  element  heat  generation 

The  data  for  each  element  must  be  placed  on  one  card.   The  nodal  point 
numbers  are  placed  in  clockwise  order  starting  with  the  number  in  the 
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upper  left  hand  corner.   Elements  may  be  generated  by  the  program  but 
there  are  some  limitations.   Only  elements  of  the  same  material  type  can 
be  generated  at  one  time.   The  heat  generation  must  be  the  same  for  the 
elements  generated  together.   Elements  must  be  generated  in  rows  or  col- 
umns,  To  generate  elements  simply  leave  out  those  element  cards  that 
you  wish  generated. 
G.   BOUNDARY  CONDITIONS 

There  are  two  types  of  boundary  conditions  that  may  be  read  into  the 
program,  convection  and  temperature  or  flux.   (An  adiabatic  boundary  con- 
dition is  a  flux  boundary  condition  with  the  flux  equal  to  zero.)   The 
following  information  is  read  into  the  program  for  convection  boundary 
conditions. 

IB         nodal  point  number 

JB        nodal  point  number 

HB         convective  heat  transfer  coefficient 

TEMPB      ambient  fluid  temperature 
The  data  for  each  portion  of  the  convection  boundary  is  put  one  one  card. 
If  time  varying  ambient  temperature  is  to  be  generated  by  the  program  the 
control  parameter  KEY  must  be  some  number  other  than  zero,  say  one.   If 
the  time  varying  ambient  temperature  is  not  generated  by  the  program  new 
data  for  each  segment  of  the  convective  boundary  must  be  read  into  the 
program  with  each  new  time  step.   In  this  case  KEY  is  zero. 

The  following  information  is  read  into  the  program  for  temperature 
and  flux  boundary  conditions. 

NN         nodal  point  number 

TQ        temperature  or  flux 
The  data  for  the  temperature  or  flux  boundary  condition  at  each  nodal 
point  is  put  on  one  card.   If  time  varying  temperature  or  flux  is  to  be 
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generated  by  the  program  the  control  parameter "JUMP  must  be  some  number 
other  than  zero.   If  time  varying  temperature  or  flux  is  not  generated, 
by  the  program  new  data  for  each  boundary  nodal  point  must  be  read  into 
the  program  with  each  new  time  step.   In  this  case  JIMP  is  zero.   Both 
temperature  and  flux  boundary  conditions  must  be  introduced  into  the  pro- 
gram in  the  same  manner.   For  both  convection  and  temperature  or  flux 
boundary  conditions,  if  time  variations  are  used  the  logical  variable 
MTAG  must  be  true.   Also,  regardless  of  the  position  of  the  last  nodal 
point,  a  boundary  condition  card  must  be  read  into  the  program  for  it. 
If  the  last  nodal  point  is  not  a  boundary  point  then  the  data  card  con- 
sists of  the  nodal  point  number  and  a  value  of  zero  in  the  temperature 
or  flux  field. 
H.   UNITS 

The  units  used  for  each  problem  must  he  consistant  with  the  input 
data  or  the  labels  on  the  output  will  be  incorrect.   Presently  there  are 
six  different  units  systems  available.   Each  system  consists  of  a  length, 
mass,  time,  and  temperature.   Each  system  has  a  number  which  must  be  read 
in  as  a  control  parameter  (UNIT).   The  systems  are: 

1.  feet,  pound  mass,  hours,  °F 

2.  inches,  pound  mass,  seconds,  °F 

3.  feet,  pound  mass,  seconds,  °F 

4.  meters,  kilograms,  seconds,  °C 

5.  centimeters,  grams,  seconds,  °C 

6.  feet,  pound  mass,  minutes,  °F 

If  the  user  does  not  prefer  the  units  available  he  can  use  another 
system  in  his  input  and  change  the  labels  or  just  ignore  the  labels  al- 
together.  The  system  of  units  will  depend  greatly  upon  the  configuration 
of  the  problem  and  the  temperature  differences  involved. 
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I.   ERROR  EXITS 

There  are  three  error  exits  in  the  program.   The  first  occurs  if  a 
non-positive  number  of  nodal  points  is  read  into  the  program.   This  exit 
is  used  to  stop  the  program.   The  other  two  occur  in  the  nodal  point  and 
element  information  reading  sections.   If  these  statements  are  executed 
then  the  execution  of  the  program  is  terminated.   If  either  of  the  last 
two  exits  are  executed  a  statement  will  be  printed  to  indicate  which  one 
it  was.   Usually  termination  means  a  bad  data  card  or  that  the  user  was 
trying  to  generate  points  beyond  the  capability  of  the  program. 
J.   COORDINATES 

The  cartesian  system  of  coordinates  x  and  y  are  used  for  problems  in- 
volving plane  geometry.   When  axisymmetric  problems  are  considered  the 
x,©-, z  system  of  polar  coordinates  is  used.   The  analysis  of  axisymmetric 
problems  is  conducted  in  the  x,z  plane. 
K.   TEMPERATURE  OUTPUT 

The  temperature  output  can  be  modified  in  two  ways.   The  number  of 
time  increments  printed  out  can  be  limited  by  using  the  control  parameter 
INTER.   If  the  program  is  directed  to  solve  for  the  temperature  every 
minute  for  100  minutes  and  INTER  is  equal  to  25,  temperature  information 
will  only  be  printed  for  25,  50,  75,  and  100  minutes.   The  second  manner 
in  which  the  output  may  be  modified  is  at  the  end  of  each  time  sequence 
the  x,y  coordinates  of  each  nodal  point  and  the  temperature  are  printed 
if  the  control  parameter  KPUNCH  equals  one. 
L.   INPUT  DATA  PREPARATION 

The  following  sections  describe  the  manner  in  which  the  input  data 
is  to  be  placed  on  data  cards.   The  data  is  prepared  in  groups  and  the 
final  data  deck  is  assembled  by  putting  the  data  groups  together  in  order 
of  their  listing  below. 
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1.   Control  Information 

(a)  Title  Card  (18A4):   72  arbitrary  spaces  for  an  alphameric 
problem  label  punched  in  columns  1-72 

(b)  Units  Card  (12M)  :   Alphameric  unit  labels 
Columns     1-12    length 

13-24    mass 
25-36    time 
37-48    temperature 
Labels  should  be  left  justified. 

(c)  Control  Information  Card  (1115,  F10.0) 


Columns 


1-5 


NUMNP 


6-10    NIMEL 
11-15 

16-20    NUMMAT 


21-25 


NSEQ 


26-30    KAT 


31-35 
36-40 
41-45 
46-50 
51-55 
56-65 
"I"  formats  must  be  right  justified. 
2.  Material  Properties 

(a)  Material  Properties  Card  (5D10.3) 
Columns     1-10    XCOND     kxx 


number  of  nodal  points 
number  of  elements 
NUMCBC    number  of  convection 
boundary  conditions 
number  of  materials 
number  of  time  sequences 
type  of  geometry 
print  interval 
units  system  (integer) 
spacial  distribution 
temperature  conditions 
convection  conditions 
initial  temperature 


INTER 

UNIT 

KPUNCH 

JUMP 

KEY 

TO 
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11-20 

YCOND 

kyy 

21-30 

XYCOND 

kxy 

31-40 

SPHT 

fc 

41-50 

HX 

heat  generation 

One  card  for  each  material. 
Time  Sequencing  Information 
(a)      Time  Sequencing   Card    (2L5,I10,D10. 3) : 


Columns 


1-5 


6-10 


LTAG 


MTAG 


11-20  ITAG 

21-21  TAG 

One  card  for  each  time  sequence. 

Nodal  Point  Information 

(a)   Nodal  Point  Card  (2I5,3F10.0) : 


temperature  dependent 
properties  (true  or  false) 
time  dependent  boundary 
conditions  (true  or  false) 
number  of  time  steps 
time  increment 


Columns 

1-5 

N 

nodal  point  number 

6-10 

K0DE 

temperature  boundary 
condition 

11-20 

X 

x  coordinate 

21-30 

Y 

y  coordinate 

31-40 

T 

initial  temperature 

One  card  for  each  nodal  point  not  generated  by  the  program. 

Element  Information 

(a)   Element  Card  (6I5,D10.3): 

Columns     1-5     N         element  number 

6-10    IX(N,1)    nodal  point  number 
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11-15    IX (N, 2)    nodal  point  number 
16-20    IX(N,3)    nodal  point  number 
21-25    IX (N, 4)    nodal  point  number 
26-30    IX(N,5)    material  number 
31-40    HTGEN     element  heat  generation 
One  card  for  each  element  not  generated.   Nodal  points  are  num- 
bered clockwise  starting  in  the  upper  left  hand  corner. 

6.  Convection  Boundary  Conditions 

(a)   Convection  Boundary  Card  (215, D10. 3,F10. 0) 

Columns  1-5  IB  nodal  point  number 

6-10         JB  nodal  point  number 

11-20         HB  convective  heat   transfer 

coefficient 
21-30    TEMPB     ambient  fluid  temperature 
One  card  for  each  convection  boundary  condition. 

7.  Temperature  and  Flux  Boundary  Conditions 

(a)   Temperature  or  Flux  Boundary  Card  (I10,D10.3) 

Columns     1-10    NN        nodal  point  number 
11-20    TQ        temperature  or  flux 
One  card  for  each  temperature  or  flux  boundary  point.   A  tempera- 
ture or  flux  boundary  card  for  the  last  nodal  point  must  be  read  even 
if  it  isn't  a  boundary  point. 

8.  Multiple  Problems 

The  above  cards  comprise  a  data  deck  for  one  problem.   Multiple 
problems  may  be  solved  by  putting  their  data  decks  together.   The  only 
limitation  on  this  procedure  is  computer  time.   It  is  important  that 
there  be  three  blank  cards  after  the  last  problem  data  set  in  the  data 
deck.   This  is  to  terminate  operation  of  the  program. 
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RECOMMENDATIONS 

Since  the  number  of  problem  types  that  can  be  solved  by  this  program 
is  limited  mostly  by  the  imagination  of  the  user  it  is  highly  recommended 
that  this  engineering  tool  should  be  brought  to  the  attention  of  those 
students  studying  in  the  area  of  unsteady  heat  conduction.   Only  through 
complete  familiarization  with  the  program  can  the  user  make  use  of  its 
great  potential  to  solve  problems  involving  complex  configurations  and 
boundary  conditions. 

Recommendations  for  further  study  in  this  area  are: 

1.  Expansion  of  the  finite  element  method  to  three  dimensional  con- 
figurations. 

2.  Revision  of  the  boundary  condition  segments  of  the  program  to 
make  use  of  random  access  disks  for  temporary  storage. 

3.  Construction  of  a  supplementary  program  that  will  punch  input 
data  cards  for  complex  time  dependent  boundary  conditions. 

4.  An  isothermal  contour  plotting  subroutine  capable  of  being  used 
in  conjunction  with  the  present  NPS  Computer  Center  subroutine  DRAW. 

5.  Use  of  the  program  in  conjunction  with  the  NPS  Computer  Center 
IBM  Optical  Display  Unit,  in  order  to  observe  transient  heat  conduction 
differences  after  minor  alterations  in  mesh  configuration. 

6.  Use  of  the  program  in  conjunction  with  existing  finite  element 
programs  in  the  area  of  stress  analysis  to  obtain  thermal  stresses. 
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APPENDIX  A 

In  this  section  some  solutions  to  different  configurations  will  be 
presented  as  examples  of  possible  problem  types  that  may  be  solved  with 
this  program.   Also  the  results  of  these  examples  are  compared  with 
analytic  solutions  of  the  same  configurations.   In  order  to  avoid  con- 
fusion the  same  material  nickel  (99.9%)  and  the  same  initial  temperature, 
100°F,  are  used  in  each  problem. 

The  material  properties  of  nickel  are: 

a  =  556  lbm/ft3 
c  =  .1065  BTU/lbm-°F 
<^T  ■  .882  ft2/hr 
Example  1. 

This  is  a  problem  in  one  dimensional  heat  conduction.   It  could  also 
be  considered  a  problem  of  heat  conduction  in  a  semi-infinite  plate  by 
assuming  that  lines  parallel  to  the  direction  of  heat  flux  are  adiabatic 
boundaries.   Because  of  this  assumption  we  only  need  one  row  of  elements 
across  the  plate  which  is  one  foot  wide  (see  figure  A-l).   The  top  and 
bottom  sides  of  the  elements  are  adiabatic  boundaries.   Also  the  end 
boundary  conditions  are  adiabatic  on  the  left  and  are  constant  0°F  on  the 
right.   The  nodal  points  are  numbered  from  left  to  right,  starting  with 

the  top  row.   The  elements  are  also  numbered  from  left  to  right.   The 
computer  solution  provides  answers  at  distinct  time  intervals  while  most 
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Figure  A-l  One  Dimensional  Heat  Conduction 
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analytic  solutions  are  given  in  terms  of  the  dimensionless  time  parameter 

o 

Fourier  number,  NpQ.   The  Fourier  number,  Npo,  equals  0L- t/L  ,  where  °^  T 

is  the  thermal  diffusivity  of  the  material.   The  comparisons  will  be  made 
at  various  Fourier  numbers.   The  analytic  solution  was  obtained  from 
page  98  of  reference  (2J. 
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Cars law  and 
Jaeger 

100 
99 
96 
85 
54 

95 
93 
82 
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Example  2. 

This  example  has  the  same  mesh  configuration  and  boundary  conditions 
as  the  first  example.   The  initial  temperature  distribution  is  changed 
to  the  form  T  e  100  -  50x°F.   The  analytic  solution  was  obtained  from 
page  98  of  reference  (jfj. 
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70 
60 

99 
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Finite 
Element 

94.0 
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79.8 
68.0 
53.0 

89.0 

85.6 
76.9 
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36 

80 
76 
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50 
28 


62.1 
37.0 

80.0 
76.0 
66.3 
49.7 
27.0 


Example  3. 

This  example  is  a  solution  for  two-dimensional  heat  conduction.   The 
configuration  is  a  square  plate  with  an  initial  temperature  of  100°F  and 
boundary  conditions  of  0°F  on  all  sides.   Using  symmetry  the  user  need 
only  examine  one  fourth  of  the  plate,  i.e.  one  corner  (see  figure  A- 2). 
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Figure  A- 2  Two  Dimensional  Heat  Conduction 

The  nodal  points  are  numbered  from  left  to  right  and  from  the  bottom  row 
to  the  top.   The  elements  are  numbered  in  the  same  manner.   In  this  case 
the  problem  consists  of  a  plate  with  two  sides  at  0°F  and  the  other  two 
insulated  for  an  adiabatic  boundary  condition.   The  analytic  solution  is 
obtained  by  using  the  solution  for  one  dimensional  heat  conduction  in  a 
semi-infinite  plate  and  the  method  of  multiplicative  superposition. 
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The  analytic  solution  was  obtained  from  page  118  of  reference  |_3 J . 

Npo 


Nodal 
Point 

8 
15 
22 
29 
36 

8 
15 
22 
29 
36 

8 
15 
22 
29 
39 


01 
01 
01 
01 
01 

05 
05 
05 
05 
05 


Chapman 


27.0 

70.6 

92.0 

97.0 

100.0 

5.75 

21.2 

40.4 

54.0 

59.3 

1.96 

7.3 

14.4 

19.4 

22.1 

Finit 

:e 

Element 

30. 

0 

70. 

6 

90. 

5 

97. 

3 

98. 

8 

6. 

29 

22. 

2 

40. 

8 

55. 

1 

60. 

4 

2. 

23 

8. 

1 

15. 

2 

21. 

1 

23. 

3 

Example  4. 

In  this  problem  the  same  mesh  configuration  as  the  previous  example 
is  used.   The  boundary  conditions  are  changed  so  that  the  sides  are  adia« 
batic,  the  top  is  at  a  constant  100°F.   The  initial  temperature  is  100°F, 
By  taking  an  extremely  large  time  increment,  that  is  of  n  magnitude 
greater  than  the  maximum  number  of  significant  figures  available  in  the 
computer,  say  10  ,  the  first  time  step  will  yield  the  steady  state 
solution  to  the  problem.   This  is  made  possible  because  '.he  t  hue  depend- 
ent contributions  to  the  coefficient  matrix  and  the  forcing  vector  are 
multiplied  by  the  inverse  of  the  time  increment  which  in  this  case  is 
essentially  zero  (10*   ).   The  solution  to  this  problem  is  a  linear 
temperature  distribution  and  the  results  below  were  obf ained  hy  using  a 
time  increment  of  10  '  minutes. 
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Finite 
Element 

0.0 
20.0 


62 


15 
21 
27 
33 


40.0 

60.0 

80.0 

100.0 


63 


APPENDIX  B 


A 

ACALTR 

B 

CONDI 

CONDJ 

CONDJ 

H 

HX 

INTER 

IT  AG 

IX 

JUMP 
KAT 
KEY 
KODE 

KPUNCH 

LPRINT 

LTAG 

MB  AND 

MTAG 

MTYPE 

NSEQ 

NUMCBC 

NUMEL 

NUMMAT 


coefficient  matrix 
dummy  logical  variable 
effective  forcing  vector 


kx  in  an  element 


ky  in  an  element 


kXy  in  an  element 


convection  boundary  layer  heat  transfer  coefficient 

heat  generation  in  a  material 

output  print  interval 

number  of  time  steps 

k  =  1  -  4  nodal  point  numbers 

k  *±  5  material  number 

temperature  and  flux  boundary  generation 

geometry  type 

convection  boundary  generation 

temperature  or  flux  boundary  parameter 

spacial  temperature  distribution  parameter 

print  parameter 

logical  variable 

half  band  width 

logical  variable 

material  number 

number  of  time  sequences 

number  of  convection  boundary  conditions 

number  of  elements 

number  of  materials 
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NUMNP  number  of  nodal  points 

SPHT  the  product  £c 

T  nodal  temperature 

TPG  time  increment 

TEMP  ambient  fluid  temperature 

TIME  time  variable 

TITLE  identification  title 

TMEAN  mean  temperature  in  an  element 

TO  initial  temperature 

UN  unit  labels 

UNIT  unit  system  parameter 

X  nodal  point  x  coordinate 

XCOND  kx  in  a  material 

XMEAN  mean  x  coordinate  of  an  element 

XYCOND  kXy  in  a  material 

Y  nodal  point  y  coordinate 

YCOND  L  in  a  material 

YMEAN  mean  y  coordinate  of  an  element 
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APPENDIX   C 

* 

Title  and  Labels 

THESIS    CHECK    PRCBLEM    HEAT    CONDUCTION    IN    A    RECTANGLE 
FEET  POUND    MASS       MINUTES  FAHRENHEIT 


Control  Information 
36         25  C  1  1  C  1  6  1  0  ClOO. 0 


Material  Properties 
♦  8.66  0D-01  +  8.660D-01  +  COOOD  +  00  +  5.9  300  +  01  +  O.OOOD  +  0  0 


* 
Time  Sequencing   Information 

FALSEFALSE  250+1 . 0CCC-C1 


Nodal   Point   Information 


1 

ICO 

CO 

100.0 

2 

10.1 

0.0 

100.0 

3 

10.2 

CO 

100.0 

4 

10.3 

CO 

100.0 

5 

10.4 

0.0 

100.0 

6 

10.5 

CO 

100.0 

7 

10.0 

0.1 

100.0 

12 

00.5 

0.1 

100.0 

13 

10. 0 

0.2 

100.0 

18 

00.5 

0.2 

100.0 

19 

10.0 

0.3 

100.0 

24 

00.5 

0.3 

100.0 

25 

10.0 

0.4 

100.0 

3C 

00.5 

0.4 

100.0 

31 

10.0 

0.5 

100.0 

36 

00.5 

0.5 

100.0 

Element   Information 

1 

7           8 

2 

1 

1+0. 0000+00 

5 

11         12 

6 

5 

1+0.0000+00 

a. 

13         14 

8 

7 

1+0.0000+OC 

1C 

17         18 

12 

11 

1+0.0000+00 

11 

19         20 

14 

13 

1+0.0000+00 

15 

23         24 

18 

17 

1+0.0000+00 

16 

25         26 

20 

L9 

1+0.0000+00 

20 

29         30 

24        ; 

23 

1+0.0000+00 

21 

31         32 

26 

25 

1+C. 0000+00 

25 

35         36 

30 

?c 

1+0.0000+00 

Boundary  Conditions 

l+O.OCCD+OO 
2+COCCD  +  OO 
3+C.OOCD+OO 
4  +  COOOD+OO 

5+o.occr+oo 
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6  +  0.CCCO00 

7+c.occr>oo 

13+C.OCCD+CO 

i°*o.occr-*oo 

25+C.OCCC+OO 
3UC.CCCP  +  CC 
36+0.0CCC+00 


NOTE:   Starred  Lines  are  not  Data  Cards 
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